Dynamical precursor of nematic order in a dense fluid of hard ellipsoids of revolution. 
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We investigate hard ellipsoids of revolution in a parameter regime where no long range nematic 
CIh' order is present but already finite size domains are formed which show orientational order. Domain 

formation leads to a substantial slowing down of a collective rotational mode which separates well 
from the usual microscopic frequency regime. A dynamic coupling of this particular mode into all 
qq ' other modes provides a general mechanism which explains an excess peak in spectra of molecular 

fluids. Using molecular dynamics simulation on up to 4096 particles and on solving the molecular 
mode coupling equation we investigate dynamic properties of the peak and prove its orientational 
origin. 

c ; 

Ch . PACS numbers: 64.70.Pf, 61.30.Cz, 61.20.Gy, 61.20.Lc, 61.25Em, 61.43.Fs 

If the temperature of a simple liquid is decreased or its density increased local correlations tend to grow. This 
tendency can be experimentally detected as a growing peak of the static structure factor S q at a wave vector q — 
Qmax — 27r/Ar at; , where Ar av is an average nearest neighbor distance. The growth of the static structure factor at 
Qmax causes the dynamic, to slow down on local length scales, a phenomenon known in neutron scattering experiments 
• as de Gennes narrowingcl. Hansen and VerletEl observed in simulations, that a large class of simple liquids crystallize, 
if the peak of the static structure factor at q max becomes larger than 2.80. The growing peak may therefore be 
Q interpreted as a precursor phenomenon of the first order crystallization transition. But since the local dynamics is 
slowing down, it can happen that the crystalline minimum of the free energy is not reached during the time S' 
the experiment. The system becomes supercooled and begins to exhibit the rich dynamics of supercooled liquid: 
Dense fluids of hard spheres are a good experimental example, where this physics can be studied in great detail. 

A simple example of a molecular liquid with translational and rotational degrees of freedom is obtained by deforming 
the hard spheres towards hard ellipsoids of revolution. They can be characterized by the diameter of the minor axes 
d and the aspect ratio Xq. This is the relation between the major and minor axes of the ellipsoid. In the following 
discussion d is always set to d — 1, i.e. length will be given in units of d, and we will address only prolate systems 
with aspect ratios larger than one. However we expect qualitatively similar results for oblate systems, replacing Xq 
by 1/aJ3. If the-aspect ratios of the ellipsoids are very large, the theory of Onsager (originally for hard sphero- 
cylinders) is validEI. This theory predicts a weakly first order phase transition from an isotropic fluid to a nematic 
upon increasing the packing fraction (p = Xopn/6, with p being the number density. Such a transition can be found 
from computer simulationsE! for aspect ratios down to Xq k, 2. The weakly first order phase transition also has a 
precursor in the static structure. The autocorrelation function S^(q — > 0) of the collective fluctuations of the second 
I ■ spherical harmonics, measuring the quadrupolar order of the molecules is strongly increasing, when the transition is 
approached. This value of the correlation function is proportional to the optical Kerr constant and related to the 
correlation function of the nematic order parameter. Theoretically it is possible to supercool also the weakly first order 
phase transition of the nematic instability!^ although to our knowledge it has not been observed yet. The increase of 
^22(1 ~ > 0) when increasing the packing fraction towards the isotropic to nematic transition leads, analogous to the 
de Gennes narrowing phenomena, to a slowing down of the correspo nding dy namic orientational correlation function. 
In extreme cases such a collective mode can drive a glass transitionOt3ll3llil. 

Our molecular dynamics (MD) simulation is placed in the regime where an increase of the peaks of the static 
structure factor, and a formation of nematic ordered domains, are both present. This has been achieved by simulating 
a system of N nearly hard ellipsoids with an aspect ratio of Xq = 1.8, where N has been chosen to be 4096 and 
512. Nearly hard ellipsoids are defined by a pair potential w(l,2) = 4e(s~ 12 — s~ 6 + |) for s < 2 1 / 6 , v(l,2) = 
otherwise. Here s = (s(l,2) — er(l,2) + d)/d is a scaled and shifted separation where the 'diameter' c(l,2) is a 
standard approximation^ to the contact distance of ellipsoids 1 and 2 at given relative orientation, and r(l, 2) is the 
center-center distance. The potential is purely repulsive and varies quite strongly over a short r. ang e of interparticle 
separation; it has been used, for higher elongations Xq, in the study of nematic liquid crystalalaO. We choose an 
energy parameter e = 1, and molecular mass M = 1. In each case the system was equilibrated at a temperature 
ksT = 1, and then constant-energy molecular dynamics were carried out using timesteps in the range St =0.003-0.005; 
typical run lengths were 10 6 timesteps. 
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From the simulation the molecular density-density correlation functions S™,(q,t) are obtained. 

S%(q,t) = ^(pt m (q,t)Pl'm(q,0)) (1) 

The indices l,l',m are the indices of the spherical harmonics, the brackets (...) denote a statistical average and 
Pim(q,t) = V^rri 1 J2n=i e ^v{iffin(t)}Yi m (p, n (t)) is the time dependent tensorial density fluctuation, where x n (t) and 
f2„(i) are position and orientation of the molecule n, respectively. We also have chosen a specific coordinate system, 
the q-frame, where the z-axis points along the wave-vector q. In this particular coordinate system all correlation 
functions depend on I and I' and only on a single helicity index m and on the absolute value q = |q|. The static 
structure factors are the equal time correlation function and reads 

S%(q) = \imS%(q,t) (2) 

For all theoretical considerations shown in this work we truncate the correlation functions at l max = rn max = 2. 

The MD simulation is first used to generate the static structure factors. In Fig. [l] for a packing fraction of <p = 0.575, 
below any transition the Sq Q (q) correlator is plotted together with the S^iq) correlator. The dense liquid manifests 
itself in the low compressibility k = kq lim^o £>oq(q)- k o is the compressibility of an ideal gas. Important there is the 
well pronounced peak at q = q ma x = 6.5, discussed above. This peak is the static manifestation of the cage formed 
by the nearest neighbors which is about to trap the particles. At higher densities it drives a glass transition. In the 
^22(9) correlator shown in the same plot an additional peak occurs at q — 0. The reader should also note that the 
peak at q = 0, which is a precursor of nematic order is accompanied by a second peak at q ~ q m ax and a third peak 
at about the same position as the third peak in S® {q). The half width of the peak at q = is approximately unity 
Aq w 1 and indicates the formation of orientationally ordered domains with a diameter of D = In j Aq ps 6. Such a 
domain formation is still far away from a true nematic transition. Yet, the S^iq) correlator reaches already a value 
of 1.9 at q = 0. The static structure factors are used as input for the_equations of the molecular mode coupling theory 
for the dynamic structure factors of supercooled molecular liquidsEj. Some details of this theory and an application 
to dipolar hard spheres can be found inliia. For the comparison with the simulation it is important to remember, 
that the mode coupling theory of supercooled liquidSi^y^estimates the slowing down of structural relaxations upon 
decreasing the temperature or increasing the densityE9u2rE3'c3c3. To allow for comparison with the simulations we 
have to rescale the packing fractions for the MMCT such, that the time scale of structural relaxations of the MMCT 
corresponds to the time scale of the structural relaxation in simulations. 

The static structure factors in dense liquids determine directly the microscopic frequencies measured in experiments 
of the respective dynamic structure factors. Neglecting the coupling to heat fluctuations in the liquid, the matrix of 
microscopic frequencies is given by 



(Bif)w = (f^+Ki+i) 1 ^-) {g-\<m, (3) 

where / is the moment of inertia of the ellipsoid with respect to a rotation on one of the minor axes and M the 
mass. For I = I' = Eq. (J3J) describes the longitudinal (isothermal) phonon frequency with its linear dispersion for 
small q- values; ^00(9) = cq with the isothermal sound velocity of the liquid c. For 1,1' 7^ and q — > the static and 
dynamic structure factors become diagonal in /, /' and independent of m for short range interaction in general and 
the hard ellipsoids in particular. Therefore the microscopic frequencies in (^) reduce for 1,1' ^ and q = to 

<tt«ftS-(l(< + l>*£)j g5 ^ W *' (4) 

For the following it is important to note, that large values of the static structure factors lead to small microscopic 
frequencies and the corresponding slow modes have a strong intensity proportional to the corresponding static structure 
factor. This is the origin of the de Gennes narrowing phenomenon in the center of mass motion at wave vectors around 
q_= qmax- The increase of the static structure factor S^(q = 0) shown in Fig.|l|, will similarly yield, due to the relation 
(Q) , a slowing down of the corresponding optical I = 2 mode. 

As soon as the structural relaxation becomes much slower than the microscopic motion, the microscopic frequencies 
become renormalized since the system is in a glassy state on the time scale of the microscopic motion. For molecular 
liquids also a splitting of the microscopic frequencies occur in addition to the renormalization due to rotation trans- 



lation coupling^ ^3. But the slowing down of the rotational mode due to the occurrence of locally ordered domains 



is preserved also under renormalizationc3. 
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Nonlinear memory effects will couple all these modes. Therefore the described slow but intense modes will influence 
the spectra at different wave vectors and different rotation numbers I in the following way. The center of mass motion 
on a spatial scale of about Ar « l/q m ax will be the primary reason for a trapping of a particle in the cage of its 
nearest neighbor which leads to the bifurcation scenario described by MCTEl Indirectly the precursor of nematic 
order seen in the static correlation function S^(q) for q = can support the tendency towards glass formation via 
the mentioned induced local spatial order, which is manifested as a peak in S^iq) at g ^ q m ax- The fact that there 
is a mode of strong oscillator strength and low frequency, S^(q ~ 0, i), can lead to a dynamic signature of this mode 
in other modes. It can be recognized by varying the inertia moment /. If and only if an excitation in a density 
spectrum at some wave vector q ^ is due to the rotational mode at q — 0, its position will shift according to Eq. (Q) 
proportional to l/VJ. The moment of inertia can be manipulated independently from the mass of the ellipsoids in 
the MD simulation and the numerical solution of the MMCT equations. Even extremely high values can be generated 
which can not be realized in a real systems. 

In Fig. U the time dependent center of mass correlation function /, I' = at q = q ma x has been plotted. For the 
moment of inertia of an ellipsoid with homogeneous mass distribution Iq = (1 + Xq)/20, we have fixed two densities 
if Sim = 0.66, ipMMCT = 0.56) which lead to approximately the same time scale for the structural relaxation in the 
simulation and the numerical solution of the MMCT, respectively. Then, the moment of inertia was increased to 2.5, 
10 and 25 times Jo. The simulation and the theoretical calculation exhibit a shift to longer relaxation times, which 
are of the same order in both cases. In addition to a shoulder at longer times a plateau develops when the inertia 
moment is increased. This demonstrates clearly the considerable qualitative influence of the rotational motion on the 
center of mass motion. To study this effect and its dependence on wave vectors more quantitatively it is instructive 
to calculate the susceptibility spectrum x"oo(9' w ) = ^"oofe UJ )- 

In Fig. | the susceptibility spectrum ~x" Q0 (q : io) of the MMCT result and the simulation are shown for two different 
wave vectors. Note that there are no adjustable parameters. For both wave vectors the wave vector dependent 
position of the high frequency excitations agree very well. The splitting of the high frequency band, seen for q — 2.71 
in the MMCT calculation and mentioned above is smeared out in the simulation. This may be understood as follows. 
Since we did not want to bias our calculation with a theory for the microscopic bare damping, we assumed a free 
oscillatory microscopic motion with frequencies given by (H). The effect of smaller damping can also be seen in 
the time dependent functions (see Fig. g) as a tendency towards oscillation in the plateau regime of the MMCT 
results. By switching on a phenomenological microscopic constant damping term in the MMCT equations the 
split microscopic frequency can be forced to merge in one single excitation as in the simulation. In addition to the 
microscopic excitation an additional almost wave vector independent hump can be seen in the MD - and the MMCT 
result at about one decade below the microscopic frequency band, before the a - relaxation peak starts to grow at 
smaller frequencies. This hump is due to the additional plateau in Fig. ^|. But the positions of the peak do not agree 
very well quantitatively The reason for this is not completely clear. The effect can partially be accounted for by the 
missing microscopic damping in the MMCT calculation, since the first minimum below the microscopic band in the 
simulation is already shifted to smaller frequencies because of damping. In addition, since the memory kernels in the 
MCT - equations, which generate the coupling between translational and rotational modes and therefore the hump, 
are approximate, only, we can in general not expect a quantitative agreement of both positions. 

The variation of the peak position of the low lying excitation in the MMCT - result with changing the moment 
of inertia can be seen in Fig. |^. The susceptibility spectra x"w (qmax, at the wave vector q max — 6.5 ~ exhibit 
the anomaly between the structural relaxation peak and the microscopic excitations for all allowed combinations of 
l,l',m, either as an additional peak or an additional shoulder. Especially in the I = I' = 2 susceptibility spectra 
a well pronounced peak can be identified, which shifts to lower frequencies for increasing moment of inertia J and 
is therefore influenced by the orientational motion. The / = I' = spectrum demonstrates that this orientation 
anomaly is hybridized with the lowest of the split high frequency excitations for I = Iq and turns into a well defined 
shoulder for larger inertia moments. These features are present in the spectrum of the MMCT calculations and the 
MD - simulations for all susceptibilities Xw with ( u ' m ) = (000), (200), (220), (221) and (222). We show some of 
them in Fig. ||. Since the MMCT calculations can be extended to longer times compared to the simulations, the 
quality of the spectra are better for the MMCT. But is is still possible to identify in Fig. [| for all moments of 
inertia the position of the anomaly. To get better statistics, we average the value of the position tu op over five spectra 
(U'm) = (000), (200), (220), (221) and (222) and plot the results against 1/V7 in Fig. §. A well defined linear law 
through the origin is found. This prov6S our assertion that the orientational anomaly in the q — Qmax 

spectra is due 

to the q = anomaly in S^(q — 0) as explained above. 

The x'' '22(91 tJ ) susceptibility of the simulation shows an additional excitation compared to m = or m = 2, which 
is not present in the version of MMCT used in this paper. This excitation is strongly wave vector dependent (see 
Fig. |t]) and turns into a shoulder for the three smallest wave vectors. It therefore does not have the same origin as 
the orientation anomaly, which is almost wave vector independent. A possible explanation could be a coupling to 
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transverse phonon modes. Such a coupling to transverse phonon modes is expected from a more general version of 
the MMCTEZI. Unfortunately the allowed wave vector range of the simulation seems to be too small to test for the 
linear dispersion lot = &rq of the expected coupling to transverse phonons in the spectrum x" 22(^1 u> )- 

In conclusion we find for the spectra Xoo"(9' w ) ancl X^'ili w ) a reasonably good agreement between the MMCT for 
hard ellipsoids and simulations of this system. We can identify a dynamic anomaly which appears about one decade 
below the microscopic excitations. This anomaly is due to the_medium range nematic order which can be identified 
by the q = value of the static structure factor S^(q = 0). InEj an analogous anomaly for dipolar hard spheres due 
to medium range ferroelectric order could be identified. There it was argued that the anomaly behaves very much 
like the well known Bose-peak phenomenon. Similar arguments can be applied to the anomaly caused by the medium 
ranged nematic order, found here for hard ellipsoids. 
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FIG. 1: For an aspect ratio of Xo = 1.8, from a MD simulation of 4096 hard ellipsoids of revolution the SQ (q) component 
(solid line) and the S^il) (dashed line) components of the static structure are shown. The packing fraction was ip — 0.575 and 
is right at the glass transition from MMCT. 
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FIG. 2: For four different values of the moment of inertia Io (solid lines), 2.5 x Iq (dotted lines) , 10 x Io dashed lines) , 25 x Io 
(long dashed lines) the time dependent relaxation of (f>o (q = q ma w,t) = Sqo(<? = q m ax,t) / So (q = q ma x,t = 0) has been plotted. 
The upper curves are shifted by 0.1 vertically and result from a numerical solution of the MMCT equations at <p = 0.56 and 
the lower curves are the result of a MD simulation with 4096 particles at tp = 0.66. 
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FIG. 3: The susceptibility spectrum x"oo(Q' u> ) from simulation at packing fraction ip — 0.66 (solid line) and MMCT at packing 
fraction ip = 0.56 (dashed line) for two different wave vectors a) qi ~ 6.5 and b) <?2 — 2.71. 
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FIG. 4: Susceptibility spectra x"w {imax, w) at the wave vector q ma x — 6.5 from the solution of the MMCT at a packing 
fraction ip — 0.56. The plots for four values of the moment of inertia 7o (solid lines), 2.5 x Jo (dotted lines), 10 x Jo (dashed 
lines), 25 x Jo (long dashed lines) are shown. 
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FIG. 6: The average position of the anomaly in the five spectra x™'"(<7> w ) is plotted against l/v7. Shown are a) the results 
for the MMCT and b) the results for the simulation. 



8 




FIG. 7: Susceptibility spectra \" 22(1^) an d x" 22(1^) fr° m t ne simulation at small wave vectors 0.4 < q < 1.33 and at 
q = q-max and tp = 0.69. 



